Esercitazione 1

Analisi dell’onda pressoria arteriosa

Per la consegna dell’esercitazione si faccia riferimento a questo link. Mentre le librerie sono rimandate nell’icona Librerie.

Parametri di campionamento che verranno considerati nell’analisi:

# Campionamento regolare
fs <- 1000      # freq. di campionamento
Ts <- 1/fs      # periodo di campionamento

1 Costruzione onda pressoria

Un metodo per interpolare dei punti è mediante una funzione spline. Peraltro, questa permette di costruire una funzione “smussata”, priva di punti singolari nel dominio interpolato.

Inizialmente è stata costruita una spline cubica naturale, che però introduce un nuovo punto di minimo in corrispondenza di \(t \approx 0.75\,\mathrm{s}\). Dunque, per garantire che venga seguito il profilo dei punti interpolati (i.e.: è necessario seguire il più possibile l’interpolazione lineare tra i diversi punti), si possono seguire due approcci:

• infittire i punti da interpolare laddove c’è una deviazione dal comportamento desiderato, ed eventualmente attribuirgli maggior importanza;

• Suggerimento ChatGPT: costruire una spline monotona shape-preserving, che consente una monotonia locale al variare dei valori del dominio e quindi permette di conservare il profilo desiderato.

Come soluzione è stato scelto l’approccio proposto dall’intelligenza artificiale.

Funzione di interpolazione.

# Punti da interpolare
tind  <- c(0, 0.1, 0.2, 0.3, 0.5, 1)   # istanti "chiave"
tval  <- c(80, 120, 100, 105, 90, 80)  # valori ai tempi chiave

# Funzione spline cubica naturale
p_spline <- splinefun(x = tind, y = tval, method = "natural")

# Segnale di onda pressoria
Ta <- max(tind) - min(tind) # tempo di acquisizioni
N <- fs*Ta                  # numero di campioni

sig_art <- tibble(
  n = 0:(N-1),
  t = n/fs,
  p_art = p_spline(t)
)

p1 <- sig_art %>% ggplot(aes(x = t, y = p_art)) +
  geom_line(color = "blue") +
  labs(x = "Tempo (s)", y = "Pressione (mmHg)", title = "Spline cubica naturale")

# Spline monotona shape-preserving
p_spline <- splinefun(x = tind, y = tval, method = "monoH.FC")

# Segnale di onda pressoria
sig_art <- tibble(
  n = 0:(N-1),
  t = n/fs,
  p_art = p_spline(t)
)

p2 <- sig_art %>% ggplot(aes(x = t, y = p_art)) +
  geom_line(color = "blue") +
  labs(x = "Tempo (s)", y = "Pressione (mmHg)", title = "Spline monotona shape-preserving")

p1+p2

Simulando dieci battiti cardiaci, il segnale complessivo risulta come segue.

# Ripetizione per 10 battiti
Ta <- 10
N <- fs*Ta

sig_art <- tibble(
  n = 0:(N-1),
  t = n/fs,
  p_art = rep(sig_art$p_art, times = 10)
)

plot_ly() %>%
  add_lines(x = t, y = p_art, line = list(width = 2.5), name = "Interp. lineare") %>%
  add_lines(x = sig_art$t, y = sig_art$p_art, line = list(width = 1.5, dash="dashdot"), name = "Spline") %>%
  layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Onda pressoria")

Si mette in rilievo la presenza di una cuspide ad ogni intero di periodo \(T = 1\,\mathrm{s}\). Per sistemare il problema si potrebbero imporre delle condizioni sul valore e sulle derivate tra le ripetizioni in corrispondenza di quei punti. Ciononostante, considerando che non è immediato il codice e che il segnale analizzato sarà quello ricostruito mediante le serie di Fourier, si trascura questo aspetto.

1.1 Ricostruzione mediante serie di Fourier

Il segnale viene ricostruito mediante le serie di Fourier, in cui l’\(n\)-esimo elemento vale:

\[ p[n] = a_0 + \sum_{k = 1}^{K} a_k \cos[2\pi k f_0 \, n] + \sum_{k = 1}^{K} a_k \sin[2\pi k f_0 \, n] , \]

dove \(f_0\) è la frequenza fondamentale, \(K\) è il numero totale di componenti armoniche e \(a_k, b_k\) sono le \(k\)-esime componenti, risultato della scomposizione su basi orto-normali. Si osservi che la componente in DC non viene divisa per 2, perché è calcolata come la media del segnale, anziché il doppio della media (come da teoria della trasformata).

La \(k\)-esima componente avrà frequenza pari a \(f_k = k f_0 = \frac{k}{T_a}\), dove \(T_a = 10\,\mathrm{s}\) è il tempo di acquisizione. Perciò, in virtù del Teorema di Nyquist-Shannon, deve valere la seguente disuguaglianza, per la quale si calcola la massima componente:

\[ f_{max} = \frac{k_{max}}{T_a} \leq \frac{f_s}2 \Leftrightarrow k_{max} = \frac{f_s T_a}{2} = 5000, \] con \(f_s = 1000\,\mathrm{Hz}\) la frequenza di campionamento.

Ricostruzione segnale anche con i coseni.

a0 <- mean(sig_art$p_art)

# La componente DC è ortogonale alle basi seno e cose
# Per evitare errori numerici, quest'ultima viene sottratta al segnale
p_art_osc <- sig_art$p_art - a0

K <- 500 # massimo possibile: 5e3
f_fund <- 1/Ta

comps_cos <- map_dbl(1:K, \(k) {
  c <- cos_k(sig_art$t, f_fund, k)
  p_art_osc %ps% c/norm_ps(c)
}) %>% zapsmall()

comps_sin <- map_dbl(1:K, \(k) {
  s <- sin_k(sig_art$t, f_fund, k)
  p_art_osc %ps% s/norm_ps(s)
}) %>% zapsmall()

comps_sin[1:50]

 [1]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
 [8]  0.000000  0.000000  9.340245  0.000000  0.000000  0.000000  0.000000
[15]  0.000000  0.000000  0.000000  0.000000  0.000000  3.153407  0.000000
[22]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
[29]  0.000000  3.164646  0.000000  0.000000  0.000000  0.000000  0.000000
[36]  0.000000  0.000000  0.000000  0.000000  1.965583  0.000000  0.000000
[43]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
[50] -0.061118

Routine mediante reduce.

# Implementazione mediante reduce
sig_art <- sig_art %>% mutate(
  p_art_k = a0 + reduce(1:K, \(acc, k) {
    c <- cos_k(t, f_fund, k)
    s <- sin_k(t, f_fund, k)
    acc + comps_cos[k] * c/norm_ps(c) + comps_sin[k] * s/norm_ps(s)
  }, .init = rep(0, n()))
)

plot_ly() %>%
  add_lines(x = sig_art$t, y = sig_art$p_art, name = "Originale") %>%
  add_lines(x = sig_art$t, y = sig_art$p_art_k, name = "Ricostruito") %>%
  layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Ricostruzione del segnale")

2 Filtraggio del segnale

Al segnale viene addizionato del rumore bianco gaussiano.

sig_art <- sig_art %>% mutate(
  p_art_n = p_art_k + rnorm(length(n), 0, 2)
)

pp <- plot_ly() %>%
  add_lines(x = sig_art$t, y = sig_art$p_art_n, name = "Con rumore") %>%
  add_lines(x = sig_art$t, y = sig_art$p_art_k, name = "Nominale") %>%
  layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Segnale ricostruito con rumore")

pp

La scelta della frequenza di taglio per il filtraggio può essere determinata sulla base dello spettro del segnale: trascurando la componente DC, che è energeticamente predominante ma non informativa sulla dinamica, è possibile valutare il taglio in funzione del contributo energetico associato a ciascuna componente in frequenza.

Tale approccio è giustificato dal Teorema di Parseval, il quale afferma che l’energia del segnale nel dominio del tempo è equivalente alla somma dei quadrati dei moduli dei coefficienti della trasformata di Fourier. In forma discreta, coerente con l’utilizzo della FFT, vale infatti: \[ \int_{0}^T s^2(t) \, dt = \sum_{k=1}^N |Y_k|^2. \]

Nello specifico il contributo relativo è calcolato come \(e_k = \frac{|Y_k|^2}{\sum_k |Y_k|^2}\).

p.fft <- sig_art %>% select(-p_art, -p_art_k) %>%
  rename(p = p_art_n) %>%
  mutate(
    f = n/Ta,
    fft = fft(p),                           # FFT concessa perché numero intero di periodi (altrimenti applicare windowing)  
    mod = (2-(f==0)) * Mod(fft)/length(n),  # riscalatura per numero di campioni e spettro monolatero
    phase = Arg(fft)/pi*180
  ) %>% slice_head(n = floor(N/2))          # spettro monolatero

p.fft %>% plot_ly(x = ~f) %>%
  add_segments(xend = ~f, y = 0, yend = ~mod, line = list(color = "#F8766D", width = 1)) %>%
  add_markers(y = ~mod, marker = list( color = "#F8766D")) %>%
  layout(
    xaxis = list(title = "Frequenza (Hz)"),
    yaxis = list(title = "Ampiezza spettrale (mmHg)"),
    title = "Spettro monolatero"
  )

p1 <- p.fft %>% dplyr::filter(f <= 10) %>%               # solo frequenze influenti
  ggplot(aes(x=f, y=mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#F8766D") +
  geom_point(colour = "#F8766D") +
  labs(x = "Frequenza (Hz)", y = "Ampiezza (mmHg)")

p2 <- p.fft %>% dplyr::filter(f<=10) %>% slice(-1) %>%  # si esclude la componente DC
  mutate(en_rel = mod^2 / sum(mod^2)) %>%               # calcola il contributo energetico relativo
  ggplot(aes(x=f, y=en_rel)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = en_rel), linewidth = 0.2, colour = "#56B4E9") +
  geom_point(colour = "#56B4E9") +
  labs(x = "Frequenza (Hz)", y = "Contributo energetico")

p1/p2

Frequenza di taglio opportuna.

Osservando gli spettri, è ragionevole adottare una frequenza di taglio di \(f_c = 2\,\mathrm{Hz}\), in quanto la componente a frequenza \(f = 1 \,\mathrm{Hz}\) contribuisce per il 69\(\,\%\) all’energia totale di oscillazione; mentre le altre componenti sono inferiori del \(10\%\).

Dunque, a partire da questa considerazione, si trova la frequenza più adatta procedendo per tentativi.

Ai fini del filtraggio, è stato impiegato un filtro digitale online butterworth del terzo ordine. Notare che l’operazione di filtraggio è offline, in quanto non si effettua in contemporanea all’acquisizione. Tuttavia, il genere di filtro impiegato è formalmente online, essendo un sistema causale (i.e: processa il segnale elemento per elemento).

fc <- 2     # freq. di taglio
fn <- fs/2  # freq. di Nyquist
ord <- 3    # filtro del 3° ordine

# filtro digitale
flt_2Hz <- butter(ord, W=fc/fn, plane="z")

p.flt_2Hz <- sig_art %>% select(-p_art) %>% rename(p_k = p_art_k, p_n = p_art_n) %>% 
  mutate(p_flt = signal::filter(flt_2Hz, p_n))

p1 <- p.flt_2Hz %>%
  pivot_longer(
    cols = c(p_n, p_flt)
  ) %>% mutate(
    name = factor(name, levels = c("p_n", "p_flt"))
  ) %>%
  ggplot(aes(x = t, y=value, colour=name)) +
  geom_line() +
  scale_colour_discrete(labels = c(
    "p_n" = "Con rumore",
    "p_flt" = "Filtrato"
  )) +
  labs(x="Tempo (s)", y="Pressione (mmHg)", colour = "Segnale", title = latex2exp::TeX("Filtraggio a $f_{c} = 2\\,$Hz"))

Alla fine viene scelta una via di mezzo: frequenza di taglio \(f_c = 5\,\mathrm{Hz}\).

fc <- 5
flt_5Hz <- butter(ord, W=fc/fn, plane="z")

p.flt_5Hz <- sig_art %>% select(-p_art) %>% rename(p_k = p_art_k, p_n = p_art_n) %>% 
  mutate(p_flt = signal::filter(flt_5Hz, p_n))

p.flt_5Hz %>% plot_ly(x = ~t) %>%
  add_lines(y = ~p_n, name = "Con rumore") %>%
  add_lines(y = ~p_flt, name = "Filtrato") %>%
  layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Filtraggio con taglio a 5 Hz")

La funzione di risposta in frequenza del filtro risulta come segue:

# Filtro butterworth con fc = 5Hz
p1 <- ggbodeplot_digital(flt_5Hz, fs, xaxis = "omega")
p2 <- ggbodeplot_digital(flt_5Hz, fs, fmin = 1e-1, fmax = 1e2) + geom_vline(xintercept = fc, linetype = 2, color = "red")

p1+p2

Eliminazione del ritardo.

L’applicazione del un filtro online introduce un ritardo rispetto al segnale originale, infatti la fase della FdT non è sempre pari a zero. Per sopperire a tale ritardo si può filtrare avanti e indietro nel tempo: in R questo filtro a fase zero è realizzato mediante la funzione signal::filtfilt. Tuttavia, sarà necessario affrontare un problema di distorsione nelle code del segnale, causato da un padding automatico della funzione.

Questa volta si può parlare formalmente di filtro offline: infatti filtfilt implenenta un sistema non-causale, visto che necessita delle misure future per un dato istante (i.e: filtraggio indietro).

p.flt_5Hz <- p.flt_5Hz %>% 
  mutate(p_flt = signal::filtfilt(flt_5Hz, p_n))

p.flt_5Hz %>% plot_ly(x = ~t) %>%
  add_lines(y = ~p_n, name = "Con rumore") %>%
  add_lines(y = ~p_flt, name = "Filtrato") %>%
  layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Filtraggio con taglio a 5 Hz")

Per eliminare la distorsione semplicemente si scartano le code: questa operazione è lecita perché il segnale si ripete periodicamente e la parte centrale è sufficientemente rappresentativa.

# [ChatGPT] si possono scartare un certo numero di campioni in base alla memoria del filtro.
# Tuttavia, questo approccio non è abbastanza, perciò vengono tagliate le code a vista.
# NB: un filtro butterworth è IIR, perciò senza memoria, ciononostante `filtfilt` usa la lunghezza
# dei coefficienti per effettuare il padding
M <- 5*length(flt_5Hz$b)

p.flt_5Hz <- p.flt_5Hz %>%
  dplyr::filter(t > 0.029 & t < 9.75)
 
p.flt_5Hz %>%
  plot_ly(x = ~t) %>%
  add_lines(y = ~p_n, name = "Con rumore") %>%
  add_lines(y = ~p_flt, name = "Filtrato") %>%
  layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Filtraggio con taglio a 5 Hz")

3 Analisi dei residui

I residui sono calcolati come la differenza tra il segnale ricostruito e quello filtrato:

\[ \varepsilon(n) = p_K(n) - p_{filt}(n). \]

Questo permette di quantificare la bonta di filtraggio: cioè quanto il segnale filtrato segue quello nominale (ricostruito mediante Fourier).

I residui verranno calcolati su tutti e tre i segnali filtrati precedentemente.

3.1 Analisi statica

Analisi statica dei residui.

# Traslo il tempo a 0: fa comodo partire da lì
p.flt_5Hz <- p.flt_5Hz %>% mutate(n=n-n[1], t = t-t[1]) %>%  
  mutate(eps = p_k - p_flt)

p.flt_2Hz <- p.flt_2Hz %>% mutate(n=n-n[1], t = t-t[1]) %>% 
  mutate(eps = p_k - p_flt)

p.flt_7Hz <- p.flt_7Hz %>% mutate(n=n-n[1], t = t-t[1]) %>% 
  mutate(eps = p_k - p_flt)

bind_rows(
  p.flt_2Hz %>% mutate(case = "2 Hz"),
  p.flt_5Hz %>% mutate(case = "5 Hz"),
  p.flt_7Hz %>% mutate(case = "7 Hz")
) %>% pivot_longer(
  cols = c(p_k, p_flt)
) %>% mutate(
    name = factor(name, levels = c("p_k", "p_flt"))
  ) %>% ggplot(aes(x = t, y = value, colour = name)) +
  geom_line() +
  scale_colour_discrete(labels = c(
    "p_k" = "Ricostruito",
    "p_flt" = "Filtrato"
  )) +
  facet_wrap(~ case, ncol = 1, scales = "free_x") +
  labs(x = "Tempo (s)", y = "Pressione (mmHg)", colour = "Segnale")

bind_rows(
  p.flt_2Hz %>% mutate(caso = "2 Hz"),
  p.flt_5Hz %>% mutate(caso = "5 Hz"),
  p.flt_7Hz %>% mutate(caso = "7 Hz")
) %>% ggplot(aes(x = t, y = eps, colour = caso)) +
  geom_point(size = 0.6, alpha = 0.6) +
  labs(x = "Tempo (s)", y = "Residui (mmHg)", colour = latex2exp::TeX("$f_c"))

Dai grafici soprastanti è evidente come anticipare il taglio distorga il segnale da quello di riferimento.

3.2 Analisi in frequenza

Analisi in frequenza dei residui.

# Funzione per automatizzare il calcolo dello spettro
compute_fft <- function(df, Ta, N, name_sig, monolateral = TRUE) {
  
  df <- df %>% mutate(
    f = n/Ta,
    fft = fft({{ name_sig }}),
    mod = ifelse(monolateral, (2-(f==0)), 1) * Mod(fft) / length(n),
    phase = Arg(fft)/pi*180
  )
  
  if(monolateral)
    return(df %>% slice_head(n = floor(N/2)))
  
  df
}

# Nuovo tempo di campionamento (a causa dell'eliminazione delle code)
Ta <- max(p.flt_2Hz$t) + 1/fs
N <- length(p.flt_2Hz$n)

# Spettro dei segnali senza la componenti a f=0
p.flt_2Hz.fft <- compute_fft(p.flt_2Hz, Ta, N, eps) %>% slice(-1)
p.flt_5Hz.fft <- compute_fft(p.flt_5Hz, Ta, N, eps) %>% slice(-1)
p.flt_7Hz.fft <- compute_fft(p.flt_7Hz, Ta, N, eps) %>% slice(-1)

bind_rows(
  p.flt_2Hz.fft %>% mutate(caso = "2 Hz"),
  p.flt_5Hz.fft %>% mutate(caso = "5 Hz"),
  p.flt_7Hz.fft %>% mutate(caso = "7 Hz")
) %>% dplyr::filter(f < 50) %>% 
  ggplot(aes(x = f, y = mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#F8766D") +
  geom_point(colour = "#F8766D") +
  facet_wrap(~ caso, ncol = 1, scales = "free_x") +
  labs(
    x = "Frequenza (Hz)",
    y = "Ampiezza Residui (mmHg)",
    title = "Spettri dei segnali filtrati"
  )

tmp <- list(fc_2 = p.flt_2Hz, fc_5 = p.flt_5Hz, fc_7 = p.flt_7Hz)

stats <-  tibble(
  fc = c("2 Hz", "5 Hz", "7 Hz"),
  mean = map_dbl(tmp, ~ mean(.x$eps)),
  sd = map_dbl(tmp, ~ sd(.x$eps))
)

knitr::kable(stats, digits = 2, caption = "Statistiche dei residui")

Statistiche dei residui

fc	mean	sd
2 Hz	0.11	6.01
5 Hz	0.01	1.74
7 Hz	-0.01	0.71

Anche osservando lo spettro dei residui è chiaro come avere una frequenza di taglio troppo bassa porti ad allontanarsi dall’andamento nominale. Infatti la dispersione dei residui tende ad aumentare (si guardi al valore di deviazione standard).

Frequenza di taglio finale.

Dunque, per concludere, da un lato non è opportuno filtrare troppo tardi, poiché non verrebbe filtrato sufficientemente il rumore del segnale. Dall’altra, anticipare eccessivamente il taglio implica una distorsione dal segnale di riferimento, di cui un buon indicatore è la deviazione standard dei residui; questo fenomeno è dovuto a una soverchiante attenuazione di componenti in frequenza, magari proprie del segnale nominale stesso. Quindi, sulla base delle precedenti osservazioni, si sceglie una frequenza di taglio pari a \(f_c = 7\,\mathrm{Hz}\), ritenuta la più approssimativa del segnale ricostruito con Fourier.